/*
 * =====================================================================================
 *
 *       Filename:  semigmm.h
 *
 *    Description:  Semisupervised GMM
 *
 *        Version:  1.0
 *        Created:  2009年06月17日 11时16分40秒
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  Ying Wang (WY), ywang@nlpr.ia.ac.cn
 *        Company:  Institute of Automation, Chinese Academy of Sciences
 *
 * =====================================================================================
 */

#ifndef SEMI_GMM_H
#define SEMI_GMM_H
#include "gmm.h"
class SemiGMM //: public preGMM
{
public:
	SemiGMM(const NCmatrix<double> &udata, const NCmatrix<double> &ldata, 
			const NCvector<int> &labels, int K, double mlambda, Covtype covtype=Full)
	//	: preGMM(unlabeldata.column())
		: unlabeldata(udata)
		, labeldata(ldata)
		, label(labels)
		, KK(K)
		, Nu(unlabeldata.row())
		, Nl(labeldata.row())
		, dim(labeldata.row())
		, type(covtype)
		, lambda(mlambda)
		, means(KK,dim)
		, fracs(KK)
		, resp(Nu,KK)
		, lndets(KK)
		, covs(KK)
		, lmeans(KK,dim)
		, lcovs(KK)
		, Nli(KK)
	{
		if( label.size() != Nl )
			std::cerr<< "the labeled data and his labels are not consistent" <<std::endl;
		NCvector<bool> iden(Nl);
		NCmatrix<double> edata;
		NCvector<double> emean;
		NCmatrix<double> ecov;
		for( int k=0; k<KK; k++ )
		{	for( int i=0; i<Nl; i++ )
			{
			//labeldata[labes[i]] Todo: compute the means and covs of labeled data
				if( label[i] == k )
					iden[i] = true;
				else
					iden[i] = false;

			}
			edata = labeldata.slicemat(iden);
			emean = mean(edata);
			Nli[k] = edata.row();
		//	ecov = cov(edata,type);
			for( int j=0; j<dim; j++ )
			{
				means[k][j] = emean[j];
				lmeans[k][j] = emean[j];
			}
			covs[k] = cov(edata,type);
			lcovs[k] = covs[k];

		}

	}
private:
	double estep();
	void mstep();
private:
    	NCmatrix<double> unlabeldata; // unlabeled data
    	NCmatrix<double> labeldata; // labeled data
    	NCvector<int> label; // label of label data
	int KK;//number of mixture
	int Nu;//number of labeled data
	int Nl; //number of unlabeled data
	int dim;// dimension of data
	Covtype type;//covariance type
	double lambda;
	NCmatrix<double> means; // means of 
	NCvector<double> fracs; // probability of each component
	NCmatrix<double> resp; //  the probability of unlabeled data belong to each component 
	NCvector<double> lndets;
	std::vector<NCmatrix<double> > covs; // the covariance of gaussian.
	NCmatrix<double> lmeans;   // means of label data
	std::vector<NCmatrix<double> > lcovs;// covs of label data
	NCvector<int> Nli;         // number of samples in each class of label data

};

double SemiGMM :: estep()
{
	int k, n;
	double sum=0, tmploglike=0,maxnum,tmp;
//	NCvector<double> x, imean;
	for ( k=0; k<KK; k++ )
	{
		for( n=0; n<Nu; n++ )
		{
			resp[n][k] = loggausspdf(unlabeldata.rowvector(n),means.rowvector(k),covs[k],type);
		}
		for( n=0; n<Nl; n++ )
		{
			if(label[n] == k)
				tmploglike += loggausspdf(labeldata.rowvector(n),means.rowvector(k),covs[k],type);
		}
	}
	for( n=0; n<Nu; n++ )
	{
		for( k=0; k<KK; k++ ) if (resp[n][k] > maxnum ) maxnum = resp[n][k];
		for(sum=0.,k=0;k<KK; k++ ) sum += exp(resp[n][k] - maxnum);
		tmp = maxnum+log(sum);
		for( k=0; k<KK; k++ ) {resp[n][k] = exp(resp[n][k] - tmp);}
		tmploglike += tmp;
	}
	//for( n=0; n<Nl; n++ )

	return tmploglike;
}

void SemiGMM::mstep()
{
	int k,n,j,jj;
	double Nui,sum1,sum2;
	if(type == Full)
	{
		for ( k=0; k<KK; k++ )
		{
			for( Nui=0.,n=0; n<Nu; n++ )
			{
				Nui += resp[n][k];
			}
			fracs[k] = (lambda*Nui + Nli[k])/(lambda*Nu + Nl);

			//Todo: mixture of covariance
			for( j=0; j<dim; j++ )
			{

				for( sum1=0.,n=0; n<Nu; n++ )
					sum1 += resp[n][k]*unlabeldata[n][j];
				means[k][j] = (lmeans[k][j] + lambda*sum1/Nli[k])/(1+lambda*Nui/Nli[k]);
				for( jj=0; jj<dim; jj++ )
				{
					for( sum1=0.,n=0; n<Nu; n++ )
						sum1 += resp[n][k] * (unlabeldata[n][j]-means[k][j]) 
							* (unlabeldata[n][jj]-means[k][jj]);
			
					covs[k][j][jj] = (lcovs[k][j][jj]+ lambda*sum1/Nli[k])/(1+lambda*Nui/Nli[k]);
				}
			}
		}
	}
	else
	{
		for ( k=0; k<KK; k++ )
		{
			for(Nui=0.,n=0; n<Nu; n++ )
				Nui += resp[n][k];
			fracs[k] = (lambda*Nui + Nli[k])/(lambda*Nu+Nl);
			for( j=0;j<dim;j++ )
			{
				for(sum1=0.,sum2=0.,n=0;n<Nu;n++)
				{
					sum1 += resp[n][k]*unlabeldata[n][j];
					sum2 += resp[n][k]*SQE(unlabeldata[n][j]-means[k][j]);
				}
				means[k][j] = (lmeans[k][j] + lambda*sum1/Nli[k])/(1+lambda*Nui/Nli[k]);
				covs[k][j][j] = (lcovs[k][j][j] + lambda*sum1/Nli[k])/(1+lambda*Nui/Nli[k]);
			}
		}
	}
}
#endif
